{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "21c23eb5",
   "metadata": {},
   "source": [
    "## 初始化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e23f723d",
   "metadata": {},
   "outputs": [],
   "source": [
    "args_wea = ['year', 'month', 'day', 'hour', 'minutes', 'temperature', 'cloud', 'sun', 'uv', 'water', 'wet']\n",
    "\n",
    "pv_time_index = ['2018/1/2', '2018/1/22', '2018/2/10', '2018/3/1', '2018/3/20', '2018/4/10']\n",
    "load_time_index = ['2018/1/2', '2018/1/22', '2018/2/10', '2018/3/1', '2018/3/20', '2018/4/10',\n",
    "    '2018/5/15', '2018/6/20', '2018/10/1']\n",
    "    \n",
    "time_str1 = '20220101'\n",
    "time_str2 = '20220430'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b006ad77",
   "metadata": {},
   "source": [
    "## 数据读取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "# macOS\n",
    "# plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']\n",
    "# Windows\n",
    "plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']\n",
    "# 正常显示负号\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "\n",
    "# 光伏\n",
    "pv_data = pd.read_csv('./data/pv_2018.csv', encoding='gbk')\n",
    "pv_data = pv_data.set_index('时间')\n",
    "print(f\"原始数据{np.shape(pv_data)}\")\n",
    "\n",
    "# 数据清洗\n",
    "from data_format import df_clean\n",
    "df = df_clean(pv_data)\n",
    "print(f\"清洗后数据{np.shape(pv_data)}\")\n",
    "\n",
    "data = np.array(pv_data['实际功率'])\n",
    "data = data.reshape(len(data) // 96, 96)\n",
    "data = pd.DataFrame(data).unstack().reset_index()\n",
    "data.columns = ['X', 'Y', 'Z']\n",
    "data['X'] = pd.Categorical(data['X'])\n",
    "data['X'] = data['X'].cat.codes\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "surf = ax.plot_trisurf(data['X'], data['Y'], data['Z'], cmap=plt.cm.viridis, linewidth=0.2)\n",
    "_ = ax.set_xlabel('时刻', fontsize = 15)\n",
    "_ = ax.set_ylabel('日期', fontsize = 15)\n",
    "_ = ax.set_zlabel('功率', fontsize = 15)\n",
    "_ = fig.colorbar( surf, shrink=0.5, aspect=10)\n",
    "# _ = plt.savefig(work_dir + '/../Image/光伏功率曲线3D.svg', format='svg')\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "840c4b2b",
   "metadata": {},
   "source": [
    "## K-Means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9ea6f78",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(pv_data['实际功率'])\n",
    "data = data.reshape(len(data) // 96, 96)\n",
    "data = pd.DataFrame(data)\n",
    "\n",
    "n_cluster = range(1, 7)\n",
    "km = [KMeans(n_clusters=i, random_state=9).fit(data) for i in n_cluster]\n",
    "scores = [km[i].score(data) for i in range(len(km))]\n",
    "_ = plt.plot(n_cluster, scores)\n",
    "_ = plt.show()\n",
    "# 将分类标签列添加至数据表中\n",
    "data['cluster'] = km[4].predict(data)\n",
    "\n",
    "# 绘图\n",
    "# colors = ['purple', 'green', 'blue', 'pink', 'brown']\n",
    "# for i, row in data.iterrows():\n",
    "#     for j in range(5):\n",
    "#         if row['cluster'] == j:\n",
    "#             _ = plt.plot(range(96), row[data.columns[:-1]], colors[j])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a59c73fa",
   "metadata": {},
   "source": [
    "## IDW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "234a3c45",
   "metadata": {},
   "outputs": [],
   "source": [
    "from idw import plot_idw_basemap, plot_idw_3d\n",
    "\n",
    "_ = plot_idw_basemap()\n",
    "# _ = plot_idw_3d()\n",
    "\n",
    "# %run ./kriging.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab391d0",
   "metadata": {},
   "source": [
    "# 发电预测"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ae0c9b",
   "metadata": {},
   "source": [
    "## 特征选取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1633dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "data = pv_data\n",
    "_ = plt.figure(figsize=(16,9))\n",
    "# mine = MINE(alpha=0.6, c=15, est=\"mic_approx\")\n",
    "mms = MinMaxScaler(feature_range=(0, 1))\n",
    "for i, time in enumerate(pv_time_index):\n",
    "    # 归一化\n",
    "    # 定性分析\n",
    "    a = mms.fit_transform(data[data.index==time]['实际功率'].values.reshape(-1,1))\n",
    "    b = mms.fit_transform(data[data.index==time][data.columns[i+3]].values.reshape(-1,1))\n",
    "    _ = plt.subplot(3, 2, i+1)\n",
    "    _ = plt.plot(a, color='b')\n",
    "    _ = plt.plot(b, color='r')\n",
    "    _ = plt.xlabel('时刻')\n",
    "    _ = plt.ylabel('功率/p.u.')\n",
    "    _ = plt.title(f'实际功率与{data.columns[i + 3]}')\n",
    "_ = plt.tight_layout()\n",
    "_ = plt.show()\n",
    "\n",
    "# # 输入变量时长确定\n",
    "# # 分析时序功率和实际功率的关系，相关性系数，决定使用多少天作为预测的时间区间\n",
    "# rows = len(data)\n",
    "# for i in range(14):\n",
    "#     a = data['实际功率'][96*14:]\n",
    "#     b = data['实际功率'][96*(14-i):(rows-96*i)]\n",
    "#     print('前' + str(i) + '天MIC：' + str(metrics.normalized_mutual_info_score(a,b)) + '\\r')\n",
    "# # 结论：选5天作为时序区间\n",
    "\n",
    "# x = pv0_data.drop(['实际功率', '时刻', '月份', '时节'], axis=1)\n",
    "# x = x.rename(columns={'风向':'Humidity', '湿度':'Wind direction', '辐照度':'Irradiance',\n",
    "#     '温度':'Temperature', '风速':'Wind speed', '压强':'Pressure'})\n",
    "# y = pv0_data['实际功率']\n",
    "# x_train,x_test,y_train,y_test = train_test_split(x,y,test_size=0.2,random_state=0)\n",
    "\n",
    "# alg = xgboost(x_train, y_train)\n",
    "# _, ax = plt.subplots(1,1,figsize=(12,6))\n",
    "# _ = xgb.plot_importance(alg, title=\"\", xlabel='Degree of importance', ylabel='Influencing factor',\n",
    "#     ax=ax, grid=False)\n",
    "# _ = plt.savefig(work_dir + '/../Image/XGBoost重要程度筛选.svg', format='svg')\n",
    "# _ = plt.show()\n",
    "# # 结论：选择温度、实际辐照度、湿度、辐照度差值作为模型的输入因素"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5234ae4",
   "metadata": {},
   "source": [
    "## 训练数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cccb500",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "def get_train_data_pv(df, x_par_names, y_par_name):\n",
    "    '''\n",
    "    df: Dataframe类型\n",
    "    x_par_names: list类型, 表中用于作为x输入的列名, 例如日照强度、湿度['sun', 'wet']\n",
    "    y_par_name: list类型, 表中用于作为y输出的列名, 例如有功['p']\n",
    "    '''\n",
    "    \n",
    "    # 时间点长度\n",
    "    time_size = 5 * 96\n",
    "    # 训练数据集与验证数据集比例为4:1\n",
    "    days = int(len(df) // (96 * 5 / 4))\n",
    "    len_train = days * 96\n",
    "\n",
    "    x_train, y_train, x_valid, y_valid = [], [], [], []\n",
    "\n",
    "    # 训练集\n",
    "    train_data = df[:len_train+time_size]\n",
    "    scaled_train_data_x = np.array(train_data[x_par_names])\n",
    "    scaled_train_data_y = np.array(train_data[y_par_name])\n",
    "    for i in range(time_size, len(train_data)):\n",
    "        x_train.append(scaled_train_data_x[(i-time_size):i])\n",
    "        y_train.append(scaled_train_data_y[i])\n",
    "    x_train, y_train = np.array(x_train), np.array(y_train)\n",
    "\n",
    "    # 验证集\n",
    "    valid_data = df[(len_train-time_size):]\n",
    "    scaled_valid_data_x = np.array(valid_data[x_par_names])\n",
    "    scaled_valid_data_y = np.array(valid_data[y_par_name])\n",
    "    for i in range(time_size, len(valid_data)):\n",
    "        x_valid.append(scaled_valid_data_x[i-time_size:i])\n",
    "        y_valid.append(scaled_valid_data_y[i])\n",
    "    x_valid, y_valid = np.array(x_valid), np.array(y_valid)\n",
    "\n",
    "    print(f'x_train shape{str(np.shape(x_train))}')\n",
    "    print(f'x_valid shape{str(np.shape(x_valid))}')\n",
    "    return x_train, y_train, x_valid, y_valid\n",
    "\n",
    "\n",
    "par_names = ['辐照度', '温度', '湿度', '实际功率']\n",
    "# pv_data[par_names] = mms.fit_transform(pv0_data[par_names])\n",
    "train_df = pd.DataFrame()\n",
    "mms = MinMaxScaler(feature_range=(0, 1))\n",
    "train_df[par_names] = mms.fit_transform(pv_data[par_names])\n",
    "x_train, y_train, x_valid, y_valid = get_train_data_pv(train_df,\n",
    "    ['辐照度', '温度', '湿度'], ['实际功率'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "dd26d1f5",
   "metadata": {},
   "source": [
    "## CNN-LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51b0a9d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras\n",
    "from keras.layers import TimeDistributed\n",
    "from keras.layers.convolutional import Conv1D\n",
    "from keras.layers.convolutional import MaxPooling1D\n",
    "from keras.layers import Flatten\n",
    "from keras.models import Sequential\n",
    "from keras.layers import LSTM, Dense\n",
    "\n",
    "# reshape from [samples, timesteps] into [samples, subsequences, timesteps, features]\n",
    "n_features = 3\n",
    "n_seq = 480\n",
    "n_steps = 1\n",
    "x_train = x_train.reshape((x_train.shape[0], n_seq, n_steps, n_features))\n",
    "\n",
    "cnn_lstm = Sequential()\n",
    "cnn_lstm.add(TimeDistributed(Conv1D(filters=64, kernel_size=1, activation='relu'),\n",
    "    input_shape=(None, n_steps, n_features)))\n",
    "cnn_lstm.add(TimeDistributed(MaxPooling1D(pool_size=1)))\n",
    "cnn_lstm.add(TimeDistributed(Flatten()))\n",
    "cnn_lstm.add(LSTM(50, activation='relu'))\n",
    "cnn_lstm.add(Dense(1))\n",
    "cnn_lstm.compile(optimizer='adam', loss='mse', metrics=['mse'])\n",
    "early_stop = keras.callbacks.EarlyStopping(monitor='val_loss', patience=10)\n",
    "# fit model\n",
    "history = cnn_lstm.fit(x_train, y_train, epochs=2, verbose=1, validation_data=(x_train, y_train))\n",
    "# y_pre = cnn_lstm.predict(x_valid)\n",
    "\n",
    "# _ = mms.fit_transform(pv_data[['实际功率']])\n",
    "# y_valid_inversed = mms.inverse_transform(y_valid)\n",
    "# y_pre_inversed = mms.inverse_transform(y_pre)\n",
    "# _ = plt.plot(y_valid_inversed[:96*5], 'b')\n",
    "# _ = plt.plot(y_pre_inversed[:96*5], 'r')\n",
    "# _ = plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e21ccf00",
   "metadata": {},
   "source": [
    "## LSTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3b07c87",
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential\n",
    "from keras.layers import LSTM, Dense\n",
    "\n",
    "# LSTM\n",
    "# x_train1, y_train1, x_valid\n",
    "# return_sequences，True：输出为一个序列，默认为False，输出一个值。\n",
    "# 输入单个样本特征值的维度\n",
    "dim = x_train.shape[-1]\n",
    "# 输入的时间点长度\n",
    "length = x_train.shape[1]\n",
    "lstm = Sequential()\n",
    "lstm.add(LSTM(units=384, return_sequences=True, input_dim=dim, input_length=length))\n",
    "lstm.add(LSTM(units=31))\n",
    "lstm.add(Dense(1))\n",
    "lstm.compile(loss='mean_squared_error', optimizer='adam')\n",
    "lstm.fit(x_train, y_train, epochs=2, batch_size=96, verbose=1)\n",
    "lstm_y_pre = lstm.predict(x_valid)\n",
    "_ = mms.fit_transform(pv_data[['实际功率']])\n",
    "y_valid_inversed = mms.inverse_transform(y_valid)\n",
    "lstm_y_pre_inversed = mms.inverse_transform(lstm_y_pre)\n",
    "_ = plt.plot(y_valid_inversed[:96*5], 'b')\n",
    "_ = plt.plot(lstm_y_pre_inversed[:96*5], 'r')\n",
    "_ = plt.show()\n",
    "lstm_y_plot = lstm_y_pre_inversed[:96]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "688beb15",
   "metadata": {},
   "source": [
    "## RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4560024c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "x_train_2D = x_train.reshape(x_train.shape[0],\n",
    "  x_train.shape[1]*x_train.shape[2])\n",
    "x_valid_2D = x_valid.reshape(x_valid.shape[0],\n",
    "  x_valid.shape[1]*x_valid.shape[2])\n",
    "\n",
    "# RF随机森林\n",
    "forest = RandomForestRegressor(\n",
    "  n_estimators=10,\n",
    "  random_state=1,\n",
    "  n_jobs=-1,\n",
    "  max_features='auto',\n",
    "  min_samples_leaf=2, min_samples_split=3)\n",
    "_ = forest.fit(x_train_2D, y_train)\n",
    "forest_y_pre = forest.predict(x_valid_2D)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "da448a35",
   "metadata": {},
   "source": [
    "## ANN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf4dd522",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPRegressor\n",
    "\n",
    "x_train_2D = x_train.reshape(x_train.shape[0],\n",
    "  x_train.shape[1]*x_train.shape[2])\n",
    "x_valid_2D = x_valid.reshape(x_valid.shape[0],\n",
    "  x_valid.shape[1]*x_valid.shape[2])\n",
    "  \n",
    "# ANN\n",
    "# 多层感知机Multi-Layer Perceptron是一种特定的人工神经网络ANN\n",
    "bpnn = MLPRegressor(\n",
    "  hidden_layer_sizes=(20,),\n",
    "  batch_size=96,\n",
    "  random_state=10, learning_rate_init=0.001)\n",
    "_ = bpnn.fit(x_train_2D, y_train)\n",
    "bpnn_y_pre = bpnn.predict(x_valid_2D)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "afc522ab",
   "metadata": {},
   "source": [
    "## MLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87a28b9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "x_train_2D = x_train.reshape(x_train.shape[0],\n",
    "  x_train.shape[1]*x_train.shape[2])\n",
    "x_valid_2D = x_valid.reshape(x_valid.shape[0],\n",
    "  x_valid.shape[1]*x_valid.shape[2])\n",
    "  \n",
    "# MLR多元线性回归Multiple Linear Regression\n",
    "linear = LinearRegression()\n",
    "_ = linear.fit(x_train_2D, y_train)\n",
    "linear_y_pre = linear.predict(x_valid_2D)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "64d15487",
   "metadata": {},
   "source": [
    "## 模型对比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b9b41f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import metrics\n",
    "\n",
    "# 预测评价指标\n",
    "def get_eval_index(y, y_pre):\n",
    "    # 均方根差RMSE、平均绝对误差MAE、决定系数R2、准确率ACC\n",
    "    rmse = round(metrics.mean_squared_error(y, y_pre)**0.5, 3)\n",
    "    mae = round(metrics.mean_absolute_error(y, y_pre), 3)\n",
    "    r2 = round(metrics.r2_score(y, y_pre), 3)\n",
    "    acc = round(1 - np.mean(np.abs((y - y_pre) / 20)), 3)\n",
    "    return {'rmse':rmse, 'mae':mae, 'r2':r2, 'acc':acc}\n",
    "\n",
    "_ = mms.fit_transform(pv_data[['实际功率']])\n",
    "y_valid_inversed = mms.inverse_transform(y_valid)\n",
    "\n",
    "# ts = [6, 14, 7, 8, 23, 24]\n",
    "# ts = [i+1 for i in range(61)]\n",
    "ts = [55, 23]\n",
    "forest_y_pre_inversed = mms.inverse_transform(forest_y_pre.reshape(-1,1))\n",
    "linear_y_pre_inversed = mms.inverse_transform(linear_y_pre.reshape(-1,1))\n",
    "bpnn_y_pre_inversed = mms.inverse_transform(bpnn_y_pre.reshape(-1,1))\n",
    "# 绘制对比曲线\n",
    "for t in ts:\n",
    "    # _ = plt.figure(figsize=(16, 8))\n",
    "    print(f\"t={t}\")\n",
    "    t1 = 96*t+24\n",
    "    t2 = 96*t+80\n",
    "    temp_data = [\n",
    "        y_valid_inversed[t1:t2],\n",
    "        forest_y_pre_inversed[t1:t2],\n",
    "        linear_y_pre_inversed[t1:t2],\n",
    "        bpnn_y_pre_inversed[t1:t2],\n",
    "        lstm_y_pre_inversed[t1:t2],\n",
    "    ]\n",
    "    temp_data = np.array(temp_data).reshape(5, 56).T\n",
    "    temp_df = pd.DataFrame(temp_data, columns=['Actual power', 'EEMD-LSTM', 'RF', 'MLR', 'SVM'])\n",
    "    temp_df.to_csv(f\"./temp_data_{t}.csv\")\n",
    "    _ = plt.plot(y_valid_inversed[t1:t2], label='Actual power')\n",
    "    _ = plt.plot(forest_y_pre_inversed[t1:t2], label='EEMD-LSTM')\n",
    "    _ = plt.plot(linear_y_pre_inversed[t1:t2], label='RF')\n",
    "    _ = plt.plot(bpnn_y_pre_inversed[t1:t2], label='MLR')\n",
    "    _ = plt.plot(lstm_y_pre_inversed[t1:t2], label='SVM')\n",
    "    print(\n",
    "        f'RF: {str(get_eval_index(y_valid_inversed[t1:t2], forest_y_pre_inversed[t1:t2]))}'\n",
    "    )\n",
    "    print(\n",
    "        f'MLR: {str(get_eval_index(y_valid_inversed[t1:t2], linear_y_pre_inversed[t1:t2]))}'\n",
    "    )\n",
    "    print(\n",
    "        f'ANN: {str(get_eval_index(y_valid_inversed[t1:t2], bpnn_y_pre_inversed[t1:t2]))}'\n",
    "    )\n",
    "    print(\n",
    "        f'LSTM: {str(get_eval_index(y_valid_inversed[t1:t2], lstm_y_pre_inversed[t1:t2]))}'\n",
    "    )\n",
    "    _ = plt.legend(fontsize=12)\n",
    "    _ = plt.xlabel('Time')\n",
    "    _ = plt.ylabel('Power/KW')\n",
    "    # _ = plt.savefig(f'{work_dir}/../Image/发电预测对比曲线{str(t)}.svg', format='svg')\n",
    "    _ = plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "447525b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import metrics\n",
    "# 预测评价指标\n",
    "def get_eval_index(y, y_pre):\n",
    "    # 均方根差RMSE、平均绝对误差MAE、决定系数R2、准确率ACC\n",
    "    rmse = round(metrics.mean_squared_error(y, y_pre)**0.5, 3)\n",
    "    mae = round(metrics.mean_absolute_error(y, y_pre), 3)\n",
    "    r2 = round(metrics.r2_score(y, y_pre), 3)\n",
    "    acc = round(1 - np.mean(np.abs((y - y_pre) / 20)), 3)\n",
    "    return {'rmse':rmse, 'mae':mae, 'r2':r2, 'acc':acc}\n",
    "\n",
    "df = pd.read_csv(\"./temp_data_55.csv\")\n",
    "for col in df.columns[1:]:\n",
    "    _ = plt.plot(df[col], label=col)\n",
    "    _ = plt.legend(fontsize=12)\n",
    "    _ = plt.xlabel('Time')\n",
    "    _ = plt.ylabel('Power/KW')\n",
    "    print(col)\n",
    "    print(get_eval_index(df[col], df[\"Actual power\"]))\n",
    "_ = plt.figure(dpi=300)\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f753801e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"./temp_data_23.csv\")\n",
    "for col in df.columns[1:]:\n",
    "    _ = plt.plot(df[col], label=col)\n",
    "    _ = plt.legend(fontsize=12)\n",
    "    # _ = plt.xlim((24, 80))\n",
    "    _ = plt.xlabel('Time')\n",
    "    _ = plt.ylabel('Power/KW')\n",
    "    print(col)\n",
    "    print(get_eval_index(df[col], df[\"Actual power\"]))\n",
    "_ = plt.figure(dpi=300)\n",
    "_ = plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a28c9a0",
   "metadata": {},
   "source": [
    "## 区间预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7130c82d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_interval(arr, dis, conf_level):\n",
    "    '''\n",
    "    arr: array数组\n",
    "    dis: 分布类型, 例如'student', 'normal'\n",
    "    '''\n",
    "    interval = None\n",
    "    mean = arr.mean()\n",
    "    std = arr.std()\n",
    "    if dis == 'student':\n",
    "        interval = stats.t.interval(conf_level, len(arr)-1, mean, std)\n",
    "    elif dis == 'normal':\n",
    "        std = arr.std(ddof=1)\n",
    "        interval = stats.norm.interval(conf_level, mean, std)\n",
    "    return interval\n",
    "\n",
    "# 预测误差\n",
    "pre_error = y_valid_inversed - y_pre_inversed\n",
    "# 预测区间\n",
    "levels = [0.95, 0.90, 0.85, 0.80, 0.75, 0.70]\n",
    "intervals = []\n",
    "for level in levels:\n",
    "    interval = get_interval(pre_error, 'student', level)\n",
    "    print(f'{level}预测下界为: {interval[0]}, {level}预测上界为:{interval[1]}')\n",
    "    intervals.append(interval)\n",
    "\n",
    "ts = [6, 8]\n",
    "for t in ts:\n",
    "    down_data = [[] for i in range(6)]\n",
    "    up_data = [[] for i in range(6)]\n",
    "    t1 = 96*t+24\n",
    "    t2 = 96*t+80\n",
    "    for i in range(56):\n",
    "        down_data[0].append(y_pre_inversed[t1:t2][i][0] + intervals[0][0])\n",
    "        up_data[0].append(y_pre_inversed[t1:t2][i][0] + intervals[0][1])\n",
    "        down_data[1].append(y_pre_inversed[t1:t2][i][0] + intervals[1][0])\n",
    "        up_data[1].append(y_pre_inversed[t1:t2][i][0] + intervals[1][1])\n",
    "        down_data[2].append(y_pre_inversed[t1:t2][i][0] + intervals[2][0])\n",
    "        up_data[2].append(y_pre_inversed[t1:t2][i][0] + intervals[2][1])\n",
    "        down_data[3].append(y_pre_inversed[t1:t2][i][0] + intervals[3][0])\n",
    "        up_data[3].append(y_pre_inversed[t1:t2][i][0] + intervals[3][1])\n",
    "        down_data[4].append(y_pre_inversed[t1:t2][i][0] + intervals[4][0])\n",
    "        up_data[4].append(y_pre_inversed[t1:t2][i][0] + intervals[4][1])\n",
    "        down_data[5].append(y_pre_inversed[t1:t2][i][0] + intervals[5][0])\n",
    "        up_data[5].append(y_pre_inversed[t1:t2][i][0] + intervals[5][1])\n",
    "\n",
    "    # _ = plt.figure(figsize=(14, 8))\n",
    "    _ = plt.plot(range(56), y_valid_inversed[t1:t2], label='Actual power', color='r')\n",
    "    _ = plt.plot(range(56), y_pre_inversed[t1:t2], label='Predicted value', color='black')\n",
    "    _ = plt.fill_between(range(56), down_data[0], up_data[0], color='#00ffcc', alpha=0.3, label=f'95% confidence interval')\n",
    "    _ = plt.fill_between(range(56), down_data[1], up_data[1], color='#00cccc', alpha=0.3, label=f'90% confidence interval')\n",
    "    _ = plt.fill_between(range(56), down_data[2], up_data[2], color='#0099cc', alpha=0.3, label=f'85% confidence interval')\n",
    "    _ = plt.fill_between(range(56), down_data[3], up_data[3], color='#0066cc', alpha=0.3, label=f'80% confidence interval')\n",
    "    _ = plt.legend(fontsize=12)\n",
    "    _ = plt.savefig(work_dir + '/../Image/光伏功率区间预测' + str(t) + '.svg', format='svg')\n",
    "    _ = plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env_b",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "vscode": {
   "interpreter": {
    "hash": "a1d72e55f37f12004b9b2d82a45c015a23cdb2733c870ac1c473b30c47b0950e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
